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Junctions and interfaces consisting of unconventional superconductors provide an excellent exper¬ 
imental playground to study exotic phenomena related to the phase of the order parameter. Not 
only the complex structure of unconventional order parameters have an impact on the Josephson 
effects, but also may profoundly alter the quasi-particle excitation spectrum near a junction. Here, 
by using spectroscopic-imaging scanning tunneling microscopy, we visualize the spatial evolution of 
the local density of states (LDOS) near twin boundaries (TBs) of the nodal superconductor FeSe. 

The 7t/2 rotation of the crystallographic orientation across the TB twists the structure of the un¬ 
conventional order parameter, which may, in principle, bring about a zero-energy LDOS peak at 
the TB. The LDOS at the TB observed in our study, in contrast, does not exhibit any signature of 
a zero-energy peak and an apparent gap amplitude remains finite all the way across the TB. The 
low-energy quasiparticle excitations associated with the gap nodes are affected by the TB over a 
distance more than an order of magnitude larger than the coherence length £ a ;,. The modification of 
the low-energy states is even more prominent in the region between two neighboring TBs separated 
by a distance « 7£, a b■ In this region the spectral weight near the Fermi level (~ ±0.2 meV) due to 
the nodal quasiparticle spectrum is almost completely removed. These behaviors suggest that the 
TB induces a fully-gapped state, invoking a possible twist of the order parameter structure which 
breaks time-reversal symmetry. 


I. INTRODUCTION 

When two superconductors are in close proximity, they 
are influenced by each other via the tunneling of Cooper 
pairs. The Cooper-pair tunneling results in the flow of a 
superconducting Josephson current, which has been stud¬ 
ied for decades and is used in various superconducting 
quantum devices [I]. The Josephson current is governed 
by the phase difference of the order parameters of the two 
superconductors. Therefore, Josephson junctions con¬ 
sisting of unconventional superconductors, where the su¬ 
perconducting order parameter changes its sign depend¬ 
ing on the momentum direction, serve as a unique plat¬ 
form where novel phase-related phenomena, e.g., sponta¬ 
neous formation of half flux quanta in a tri-junction of 
cuprate superconductors Q, take place. Compared with 
the well-investigated Josephson currents, the spatial and 
energy dependence of the superconducting order param¬ 
eter and quasiparticle states around these junctions re¬ 
main to be understood. 

Recent progress in scanning tunneling microscopy 
(STM) and spectroscopy (STS) technologies opens up a 
way to directly visualize the spatial variation of the elec¬ 
tronic states in superconducting hetero-structures M- 
However, STM/STS studies on superconducting junc¬ 
tions made of unconventional superconductors are still 
demanding. There are two reasons which make it diffi¬ 
cult to study unconventional junctions. First, it is often 


challenging to artificially fabricate well-defined junctions 
of unconventional superconductors. Second, in most of 
unconventional superconductors, surfaces are not elec¬ 
tronically neutral; the resultant charge accumulation at 
the surfaces prevents STM/STS from accessing bulk su¬ 
perconducting properties. In this study, we solve these 
problems by inspecting the twin boundaries (TBs) in the 
nodal iron-based superconductor FeSe [g, 3 ]- 

The TB is a crystallographic plane in a crystal 
shared by two neighboring domains with one being 
the mirror image of the other. The TBs are often 
formed by a tetragonal-to-orthorhombic structural phase 
transition, which reduces the four-fold (C4) symme¬ 
try at high temperature to two-fold (C2) symmetry 
at low temperature. In such a case, the orthorhom¬ 
bic crystal may contain the TBs parallel to the ( 110 ) 
plane, which act as an atomically well-defined junction. 
Some unconventional-superconductor-related materials, 
such as YBa2Cii307_5, AE(Fei_ a ,Co a; )2As2 (AE: alkali- 
earth element) and NaFeAs, do form TBs upon the 
tetragonal-to-orthorhombic transition which were iden¬ 
tified by STM/STS measurements (l-doj. However, 
unavoidable surface state formation and/or insufficient 
amount of chemical doping prevent the STM/STS mea¬ 
surements to access superconductivity near TBs in these 
materials. 

FeSe (superconducting transition temperature T c « 
9 K 0) is a promising candidate for studying the ef- 
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fects of TBs on unconventional superconductivity by 
STM/STS. Among various iron-based superconductors, 
FeSe has the simplest crystal structure [Fig. 1(a)] in 
which electronically neutral two-dimensional FeSe layers 
are stacked along the c axis 0- This guarantees the per¬ 
fect cleaved surface which is electronically neutral. The 
tetragonal-to-orthorhombic structural phase transition, 
which is likely caused by the orbital ordering QM3 , oc¬ 
curs at T s ft 90 K and the TBs are spontaneously formed 
in the orthorhombic phase as illustrated in Fig. 1(b). 

Band-structure calculations show that the Fermi sur¬ 
face of FeSe consists of hole cylinders around the 
zone center and c omp ensating electron cylinders around 
the zone corner [18], |l9[. Several measurements, in¬ 
cluding penetration depth, quasiparticle interference, 
thermoelectric response Q, quantum oscillations |T~il . 
l20l . l2lj| . and angle-resolved photoemission spectroscopy 
(ARPES) I I2l I IT 12 2l reveal that the Fermi surface in the 
orthorhombic phase consists of one hole and one (or two 
hi mi) electron bands, both of which have very low 
carrier densities. The tunneling spectrum Q, tempera¬ 
ture dependence of the penetration depth down to 80 ml( 
and the residual thermal conductivity at T —> 0 [7j, all 
provide strong evidence that FeSe is an unconventional 
superconductor with line nodes in the superconducting 
gap. 

The TBs in FeSe have been studied by low-temperature 
(4.2 K) STM/STS in the films grown by molecular beam 
epitaxy and the suppression of superconductivity by the 
TBs has been reported (2[|- We performed STM/STS 
measurements at much lower temperature (ft 0.4 K) in 
vapor-grown single crystals to examine the details of the 
superconducting gap and quasiparticle excitations near 
the TBs. 


II. EXPERIMENTAL METHOD 

STM/STS experiments were conducted in a constant- 
current mode with a commercial ultra-high vacuum very- 
low temperature STM (UNISOKU, USM-1300) modi¬ 
fied by ourselves [24|]. The samples used in this study 
were high-quality bulk single crystals grown using the 
vapor transport method [26[. Superconducting tran¬ 
sition temperature defined at zero resistance is about 
9 K. These crystals are undoped and stoichiometric, en¬ 
abling us to investigate uniform and clean TBs. Sam¬ 
ples were cleaved in-situ at liquid N 2 temperature to 
prepare clean and flat (001) surfaces. Immediately af¬ 
ter cleaving, the samples were transferred to the STM 
unit kept below 10 K. We used electrochemically-etched 
polycrystalline tungsten wires for the scanning tips which 
were cleaned and sharpened in-situ by field evaporation 
using field-ion microscopy. The tunneling conductance 
g(r,E) = dIt/dV s (r,E ) reflecting the local density of 
states (LDOS) at a position r and energy E, was ac¬ 
quired by standard lock-in technique. Here, It and V s 
denote the tunneling current and the sample-bias volt¬ 


age, respectively. 

III. RESULTS AND DISCUSSION 
A. Imaging the twin boundary in FeSe 

Figure 1(c) depicts an STM image of the cleaved sur¬ 
face of an FeSe single crystal at temperature T = 1.5 K. 
The image demonstrates the extremely small concentra¬ 
tion of defects, i.e. about one defect per 5000 Fe atoms 
in the (001) plane. There is a shallow “groove” running 
along the [110] direction of the Fe lattice across which 
the unidirectional feature around the point defect is ro¬ 
tated by 7t/2, indicating that the “groove” represents the 
TB. We also observed that the elongated vortex cores @, 
which were imaged by mapping g{r , E = 0) in a magnetic 
field, are rotated by 7t/2 across the TB [Fig. 1(d)]. What 
is intriguing is that the vortices trapped at the TB are 
not elongated along the TB, demonstrating that the crit¬ 
ical current density across the TB is comparable to that 
of the bulk. The STM image of the TB at a lower bias 
voltage is not a “groove” but a “ridge” [Figs. 1(e)]. This 
suggests that the apparent corrugations near the TB are 
primarily associated with the electronic-state variations; 
the actual surface topography near the TB may be essen¬ 
tially flat. A magnified STM image near the TB is shown 
in Fig. 1(f). A regular square lattice of the top-most Se 
atoms is well maintained even in the close vicinity of the 
TB. These observations indicate that the TB in FeSe is 
an atomically sharp superconducting junction with mini¬ 
mal strain to the lattice. (A detailed argument regarding 
the absence of strain is given in Appendix A.) 

B. Local density of states across the twin 
boundaries 

We examined the LDOS evolution across the TB by 
taking a series of g(r,E) along the line indicated in 
Fig. 2(a). Here and in the following, we are interested 
only in the evolution of g(r,E ) along the x axis run¬ 
ning perpendicular to the TB leaving the y coordinate 
constant, hence g[x,E). Figure 2(b) shows an intensity 
plot of g(r,E). Individual spectra taken at representa¬ 
tive points are depicted in Fig. 2(c). At the position far 
away from the TB (I), g(r,E) exhibits a superconduct¬ 
ing gap with clear quasiparticle peaks at ft ±2.5 meV. In 
addition to this main feature, there is a shoulder outside 
of the main peaks (ft ±3.5 meV), which may represent 
multiple superconducting gaps 0 • In contrast to the case 
of fully-gapped superconductors in which g(r , E) = 0 in 
an extended E region near E = 0, g(r, E) in FeSe ap¬ 
proaches zero only for E —> 0 and apparently V-shaped, 
indicating the presence of line nodes Even right at 
the TB (III), the residual LDOS at E = 0 is negligi¬ 
bly small, indicating that the TB hardly gives rise to a 
pair breaking effect. In the vicinity of the TB, the quasi- 
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FIG. 1: (a) Crystal structure of FeSe visualized using the VESTA program [25], (b) Schematic top view of the atomic 

arrangement near the TB of FeSe (not in scale). Green filled circle and orange open circle denote top-most Se and Fe atoms, 
respectively. Se atoms beneath the Fe layer are not shown, (c) A constant-current STM image of the cleaved (001) surface 
of FeSe at 1.5 K showing the TB running from bottom left to top right. Crystallographic axes parallel to the Fe-Fe direction 
are shown by white arrows (b > a). The two insets show a magnified image of the defect (8.8 nm x 8.8 nm) in the upper-left 
or lower-right domain. Note that the pattern is rotated by n/2 between the two domains. The set-up conditions for imaging 
were 14 = +95 mV and It = 10 pA. (d) Zero-bias conductance image g(r,E = 0) at 1.5 K showing vortices. A magnetic field 
of 1 T was applied along the c axis. The tip was stabilized at 14 = +10 mV and It = 100 pA. A bias modulation amplitude 
14nod = 0.21 mV r ms was used for spectroscopy, (e) A low-bias STM image at 1.5 K taken with 14 = +20 mV and It = 10 pA. 
The field of view for (c)-(e) is the same, (f) An atomic-resolution STM image near the TB which is running vertically in the 
center of the field of view. 14 = +95 mV and It = 100 pA. 


particle peak and the shoulder associated with the su¬ 
perconducting gap diminish, and instead, sharp particle- 
hole symmetric peaks appear at E r; ±1.5 meV. In the 
crossover region (II), the 1.5 meV peak coexists with the 

2.5 meV peak, meaning that the former is not associated 
with the suppressed superconducting gap. The 1.5 meV 
peak diminishes within a distance of about 5 nm from the 
TB, which is close to the “averaged” in-plane coherence 
length £ a b r; 5 nm obtained from the upper critical field 
H c 2 (11 c) r;15 T H3- These results suggest that the 

1.5 meV peak represents the bound state induced by the 
TB. 

Another interesting observation is that low-energy 
quasiparticle excitations are modified over a very long 
distance from the TB. High-resolution g(r,E ) spectra at 
the positions of (I), (II) and (III) are plotted in Fig. 2(d). 
While the overall V-shaped behavior is maintained, the 
exact shape near the bottom of the gap depends on the 
position. In order to examine this behavior, we fit an 
empirical power-law g(r,E) oc \E\ a to the low-energy 
(l-El < 0.5 meV) spectra and plot the exponent a as a 
function of the distance from the TB at x = 0 [Fig. 2(e)]. 
Except close to the TB (]x| < t; a b) where the 1.5 meV 


peaks dominate, a increases gradually with decreasing x 
by about r: 40%. This implies the suppression of the low- 
energy quasiparticle excitations, most probably due to 
the opening of a small gap induced by the TB. The salient 
feature is that a continues to change even at |x| > 10£ o b 
(r; 50 nm), indicating an unexpectedly long-distance in¬ 
fluence of the TB. 

The long-distance TB effect on the LDOS can be seen 
in a more dramatic way in two junctions in series formed 
by two TBs. As shown in Fig. 3(a) we find an area where 
two TBs are running parallel to each other. The dis¬ 
tance between the TBs is 34 nm, which is about 7 times 
larger than £ a {,. Figure 3(b) shows the spatial evolution 
of g{r , E) across the double TBs. Individual spectra at 
representative points are plotted in Fig. 3(c). The overall 
spectral features, the 2.5 meV peak, the 3.5 meV shoul¬ 
der and the 1.5 meV peak observed near a single TB are 
all reproduced (positions I, II, and III). However, the low- 
energy spectrum taken inside the central domain (posi¬ 
tion IV) shows a striking anomaly which is absent in the 
case of a single TB. Figure 3(d) depicts g(r,E) spectra 
at low energies. It is clear that, in between the double 
TBs, there is a finite energy range where g(r, E) is almost 
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FIG. 2: (a) A constant-current STM image near a TB taken at V„ = +95 mV and I t = 100 pA. (b) Intensity plot of g(r, E ) 

along the yellow broken line in (a). V s = +20 mV, It = 100 pA and V mo d = 0.05 mV rm s. (c) Tunneling spectra at the 
representative points indicated in (a). Positions (II) and (II’) are symmetric about the TB. V s = +20 mV, It = 100 pA and 
Vmod = 0.05 mV rms . (d) High-resolution tunneling spectra at low E taken at the same positions as for (c). V 3 = +10 mV, 
It = 100 pA and Vm 0 d = 0.025 mVrms. Open symbols and solid lines denote experimental data and fitted results, respectively. 
Spectra shown in (c) and (d) are shifted vertically for clarity, (e) The exponent a ( g(r,E ) oc \E\ a ) determined from the fit to 
the experimental data in the range of |i?| < 0.5 meV plotted as a function of the distance from the TB. 


completely zero. The noticeable difference of the gap 
structure between inside and outside the central domain 
is clearly seen in Fig. 3(e), which shows the exponent a 
plotted as a function of the distance from one of the TBs; 
a is strongly enhanced in the central domain and peaks 
at the middle of the domain. The large power a w 4, 
which is ~ 3 times larger than the values at large x , 
essentially indistinguishable from an exponential energy 
dependence. This apparent large power again corrobo¬ 
rates the finite gap opening in the excitation spectrum of 
quasiparticle. 


C. Possible time-reversal-symmetry-broken state 
near the twin boundary 

The above observations, the TB-induced bound states 
at finite energies and the suppression of the low-energy 
quasiparticle excitations over a length scale much longer 
than £ a t>, suggest a novel role of the TB in an unconven¬ 
tional superconductor. Before discussing the origin of 
these anomalies, we briefly review what can be expected 
at a TB of FeSe. Recent high-resolution laser-ARPES 
measurements of FeSe indicate that the hole cylinder is 
fully gapped {27} , implying that the line nodes are present 


on the electron cylinder, as has been also inferred from 
vortex imaging |6|. Given this information, we consider 
two possible phase structures for symmetry of the su¬ 
perconducting gap across a TB as illustrated in Fig. 4, 
where either the global phase of the superconducting gap 
is fixed to the crystallographic axis [Fig. 4(a)] or is flipped 
across the TB [Fig. 4(b)]. It should be noted that the sign 
of either the nodal gap or the nodeless gap is reversed 
between the two domains in Fig. 4(a) or Fig. 4(b), re¬ 
spectively. This means that the amplitude of at least 
one of the gaps vanishes at the TB, giving rise to the 
zero-energy quasiparticle state that should appear as a 
zero-energy peak in g(r,E). This argument applies not 
only for the particular phase structure shown in Fig. 4 
but also for a general case in which nodal and nodeless 
gaps reside on multiple Fermi surfaces. 

The observed bound-state peak at 1.5 meV apparently 
contradicts this conjecture and suggests instead that the 
TB induces an additional gap component which shifts 
the position of a zero-energy peak to a finite energy. We 
point out that, as long as the induced gap is real, a sum 
of the bulk gap and the TB-induced gap reverses its sign 
at a finite distance from the TB and still gives rise to 
a zero-energy peak. However, as shown in Fig. 2(b), we 
did not observe a zero-energy peak in g(r , E) over more 
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FIG. 3: (a) Constant-current STM image with double TBs taken at Id, = +95 mV and I t = 10 pA. (b) Intensity plot of 

g(r,E) along the yellow broken line in (a). Positions of the TBs are indicated by broken lines. A low-conductance position 
at «-47 nm is due to a point defect nearby. V s = +20 mV, I t = 100 pA and V mo d = 0.05 mV rms . (c) Tunneling spectra at 
the representative points I to IV indicated in (a). V a = +20 mV, It = 100 pA and Vmod = 0.05 mV rms . (d) High-resolution 
tunneling spectra at low E taken at the same positions as for (c). V s = +10 mV, It = 100 pA and V m od = 0.025 mVrms- 
Symbols and solid lines denote experimental data and fitted results, respectively. Spectra shown in (c) and (d) are shifted 
vertically for clarity, (e) The exponent, ( g(r,E ) oc |-E|“) determined by the fitting in the range of |A| < 0.5 meV plotted as a 
function of the distance from one of the TB. The data of the single TB is shown by a red line for reference. 


than 100 nm from the TB. Thus, we speculate that the 
induced gap has an imaginary component, which means 
that time reversal symmetry is locally broken near the 
TB. In such a case, bound state peaks are located at 
finite energies E = ±A cos(<5<^/2) because the phase shift 
8ip on the TB is reduced from 7 r [28|, [2{|. Here, A is the 
amplitude of the superconducting gap. The possibility 
of the TB-induced time-reversal-symmetry-broken state 
has been argued theoretically in d -wave YBa 2 Cu 3 07_5 
with a small s-wave component pTo j] , but the experimental 
observation is still lacking. 

In order to substantiate the relevance of this scenario, 
we have calculated the spatial evolution of the LDOS 
for a model order parameter with broken time-reversal 
symmetry near TBs. The ^-symmetric order parameter 
is represented by a sum of the isotropic component Ai so 
and the four-fold nodal component A 4 ^ sin( 2 </>), 

A(x) = A iso + A 4 ^(x) sin( 2 </>), ( 1 ) 

where </> is the azimuthal angle in the momentum space; 
see Fig. 4. We assume that the global phase of the order 
parameter is fixed to the crystallographic axis, that is, 
the nodal component changes its sign across a TB as 


shown in Fig. 4(a). The spatial variation of A 4 ^ around 
a TB at x = xq is modeled by the form 

A 4 <? j(x) = A 40 lk {tanh[(x — xo)/£] cosd(x) + isin0(x)}, 

( 2 ) 

where the x axis is taken to be perpendicular to the TB, 
and A^ lk is the amplitude of A 4 ^ in the bulk. The phase 
<p of A 4 ^(x) equals 9(x) for x — Xo £ and n — 6(x) for 
— (x — xo) £. The phase variable 9{x) is assumed to 
take a nonvanishing value near the TB and exponentially 
decay with another length scale £. It is important to 
note that the characteristic length £ for the local time- 
reversal symmetry breaking can be much longer than the 
coherence length £ [30]. (The derivation of the length 
£ is given in Appendix B.) To account for low-energy 
excitations near the nodes observed in the LDOS, we 
focus on the electron cylinder with nodal gaps by setting 
the parameters Aj so = 0.2A 0 and A^ lk = 0.8A o . 

As a model order parameter with a TB at x = Xo = 0, 
we take 9(x) = ( 7 r/ 6 )sech(x/£) with £ = 5£, which gives 
A 4 ^(x = 0) = (i/2)A^ lk . The order parameter A 4 ^(x) 
is shown in Fig. 5(a), where | A 4 ^| changes with the length 
scale £ while Im(A 4 0 ) decays with the longer length scale 
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FIG. 4: Schematic illustration of the phases of the superconducting gaps across the TB shown by the red line. Top panel 
represents the iron lattice near the TB together with the momentum-space phase structure of the superconducting gaps opening 
at multiple Fermi cylinders, a hole cylinder at the center and electron cylinders at the corner of the Brillouin zone (black broken 
square). Different colors (red and blue) denote different signs of the phase. We assume that the gap node exists on the electron 
cylinder and the sign reversal is between the main lobe of the gap on the electron cylinder and the gap on the hole cylinders 
but the argument given in the text applies not only for this particular case but also for other cases. There are two possibilities: 
either the phase structure is fixed to the lattice (a) or is flipped across the TB (b). In the former case, the nodal component 
A 4 ^ should change its sign across the TB, whereas the sign of the isotropic component Ai so (either due to the fully gapped 
Fermi cylinder or associated with the C 2 symmetry of the nodal gap) should be reversed in the latter case. 


£. The phase <p abruptly changes near the TB and grad¬ 
ually approaches 0 or n. Using this order parameter, 
we calculate the spatial dependence of the LDOS within 
the quasi-classical approximation j3lj]. Figure 5(b) shows 
the global peak structure of the LDOS at representative 
points calculated with energy smearing of rj = 0.03Ao. 
Far from the TB, namely in the bulk (I), the LDOS 
has peaks at |£’| = A^ lk ± Ai so . On the TB (III), 
the peaks observed in the bulk are suppressed, and al¬ 
ternative peaks appear at E « ±0.4Ao, which corre¬ 
spond to the bound states whose energies are shifted from 
E = 0 due to the local time-reversal symmetry breaking 
in A 4 ? j. The bound-state peaks disappear at x = 3£ (II) 
since their wave functions decay into the bulk with the 


length scale £. These features of the calculated LDOS 
arising from the bound states are consistent with the 
LDOS peaks observed at Em ±1.5 meV by STM. We 
show in Fig. 5(c) the LDOS at lower energy scale which 
has been calculated with a much smaller smearing factor 
rj = 0.001A 0 . The clear V-shaped LDOS in the bulk (I) 
changes to the U-shaped LDOS upon approaching the 
TB, in agreement with the increase of the exponent a 
evaluated from the experimental data [Fig. 2(e)]. The 
low-energy LDOS is finite at x = 0 (III) and x = 3£ (II) 
because low-energy quasiparticles with momenta along 
the nodal directions of the bulk gap can linger over long 
distance and reach the TB, even though the local gap, 


| A (a:) | = [A iso + Re(A 4 ^) sin(2 cj))} 2 + [Im(A 4<? j) sin(2</>)] 2 , (3) 

does not vanish near the TB where Im(A 4 ^) ^ 0. 

We also calculate the LDOS for double TBs located at x = ±x 0 = ±3.5£, taking the model order parameter of the 
form 
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FIG. 5: (a) A model order parameter A^a;) with a TB located at x = 0. LDOS’s in the bulk (I), at x = (II), and on a TB 

(III), which are calculated with an energy smearing of g = 0.03Ao (b) and g = O.OOlAo (c). The lines (I) and (II) have offsets 
go and 0.5(?o in (b) and 0.3go and 0.15<?o in (c), respectively, where go is the density of states in the normal state at the Fermi 
energy, (d) A model order parameter A^(x) with double TBs located at * = ±3.5£. The LDOS in the bulk (I), at x = 7£ (II), 
on a TB (III), and at the middle point between double TBs (IV), which are calculated with g = 0.03Ao (e) and g = O.OOlAo 
(f). The lines (I), (II), and (III) have offsets go, ( 2/3)go , and (l/3)go in (e) and 0.3go, 0.2go, and O.lpo, in (f) respectively. 


A 40 (x) = A40 lk {tanh[(a; — £o)/£] tanh[(a; + x 0 )/£] cos 0 (a;) + *sin 0 (a;)}, ( 4 ) 


shown in Fig. 5(d). We assume that the distance 2xq 
between the TBs is in the range £ -C xo < £. The phase 
9{x) is an even function of x and takes a maximum value 
at x = 0. The global peak structure of the LDOS and 
its low-energy blowup are shown for representative points 
along the x direction in Fig. 5(e) and 5(f), respectively. 
The large peaks at \E\ « 0.4Ao on a TB (III) and the 
small peaks at \E\ ss 0.2Ao at the middle point x = 0 be¬ 
tween the TBs (IV) in Fig. 5(e) originate from the same 
dispersive mode of bound states at a TB. The calculated 
LDOS spectrum at x = 0 between the double TBs (IV) in 
Fig. 5(f) exhibits a clear energy gap extending over the 
region |i?| < 0.1A 0 , reflecting the existence of a larger 
local gap at x = 0, where the bulk low-energy quasi¬ 
particles cannot reach. We conclude that the local gap 
enhanced by the local time-reversal symmetry breaking 
near TBs over the length scale £ can explain the strong 
suppression of the LDOS between the two TBs observed 
in our STM/STS experiments. 


IV. SUMMARY 

We have reported on the visualization of the atomic 
scale variation of the quasiparticle states of the nodal 
superconductor FeSe near TBs that enforce a sign in¬ 
version at least one of the superconducting gaps open¬ 
ing on multiple Fermi cylinders. In contrast to the ex¬ 
pectation that the sign inversion generates a zero-energy 
quasiparticle bound state near the TB, the TB-induced 
quasiparticle states are not at zero but at finite ener¬ 
gies E ~ ±1.5 meV. Moreover, the low-energy excitation 
spectrum is affected by the TB over an extremely long 
distance, which is a few tens of times larger than £ a b. 
An even more dramatic change in the low-energy spec¬ 
trum has been detected in the region between double 
TBs separated by a distance ~ 7£ a &, where the quasipar¬ 
ticle weight near the Fermi energy is almost completely 
removed in the energy range \E\ < 0.2 meV. These ob¬ 
servations are qualitatively reproduced by a phenomeno¬ 
logical model which assumes that the TB induces locally 
a superconducting state that breaks time-reversal sym¬ 
metry. 























Our results suggest several important directions for 
future studies. A microscopic mechanism that induces 
the time-reversal symmetry broken state is elusive and 
should be investigated theoretically. It is also interesting 
to go beyond the quasi-classical approximation because 
FeSe is a unique superconductor whose Fermi energy is of 
the same order as the superconducting gap, placing this 
system in the BCS-BEC crossover regime [7|. Exper¬ 
iments that directly probe the time-reversal symmetry 
breaking, such as muon-spin rotation and local magne- 
tometry, are highly desirable and would give us further 
insights into the unconventional superconducting junc¬ 
tions. We anticipate that the TBs in FeSe will stimulate 
further research on the role of the phase of the supercon¬ 
ducting order parameter near the interface, which has 
been difficult to access experimentally. 
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Appendix A: Absence of lattice distortion induced 
by the twin boundary 

Although STM has a high spatial resolution, possi¬ 
ble creep in the piezoelectric scanner and/or the thermal 
drift make it difficult to estimate the small distortions 
in the topographic image. Here we utilize the so-called 
Lawler-Fujita algorithm [32| to deduce the lattice distor¬ 
tion and show that the TB-induced strain is negligibly 
small. 

First we briefly explain the principle of the methodol¬ 
ogy. The observed STM topography T(r), which mainly 
represents the top-most Se lattice, can be expressed as 


T(r) = 

T 0 [cos {q x ■ (r - u(r))} + cos {q y ■ (r - u(r))}] H-. 

(Al) 

Here, Tq is the amplitude of the Se-lattice modulation, 
q x and q y are wave vectors for the Se lattice, and • • • 
represent all other modulations. The distortions from the 
perfect lattice is described by the displacement field u(r) 
that can be regarded as a spatially varying phase of the 
q x and q y modulations. This approximation is justified 
as long as the length scale of distortions is much longer 
than the Se-Se distance as e - Standard phase-sensitive 
detection scheme can be used to evaluate u(r). By mul¬ 
tiplying T(r) and the reference signal cos (q x ■ r), we get 


T(r) cos {q x ■ r) 


n 

2 


cos (q x ■ u(r)) 

+ cos (2 q x ■ r — q x ■ u(r)) 

+ cos {{q x + q y ) ■ r — q y ■ u(r )) 


+ cos {(-q x + q y ) r-qy u(r)) 


(A2) 


All terms except the first exhibit periodic spatial modula¬ 
tions, which can be removed by low-pass Fourier filtering 
LPF {•••}. 


LPF {T(r) cos (q x •»")} = -y cos (q x ■ u(r)). (A3) 

By using the quadrature reference sin (q x • r), we get 


LPF {T(r) sin {q x ■ r)} = y sin (q x ■ u(r)) . (A4) 

Therefore, we obtain u x (r), the x component of u(r) as 


U (r) = £Se t -1 LPF |T(r) sin (q x ■ r)} 
x 27r LPF {T(r) cos (q x ■ r)} 

The y component u y (r) can also be deduced as 


u y {r) = T^tan 


-l 


LPF {T(r) sin (q y ■ r)} 
LPF {T(r) cos ( q y ■ r)} ’ 


(A6) 


A schematic model of atomic arrangement near the TB 
is shown in Fig. 6(a). We expect that the orthorhombic 
distortion affects the atomic arrangement along the y di¬ 
rection across the TB, while the periodicity along the 
x direction remains intact. In order to verify this model 
and to check if there is an additional lattice distortion, we 
calculated u x (r ) and u y (r) of the high-resolution STM 
image containing a TB running along the y direction 
[Fig 6(b)]. Reference wave vectors q x and q y were ob¬ 
tained by the Fourier analysis in the left domain. For low- 
pass Fourier filtering, we picked up only long-wavelength 
components by using a Gaussian mask with half width at 
the half maximum of 0.21(27r/as e ). Since there is a trans¬ 
lational symmetry along the TB, we average u x (r) and 
u y {r) along the y direction, yielding u x vg (x) and u y vg (x ), 
respectively. This significantly enhances the signal-to- 
noise ratio. 

Figures 6(c) and (d) show u| vg ( x) and its x deriva¬ 
tive. There is no noticeable anomaly in both u^ s (x) and 
du x vg (x)/dx, except for the smooth background associ¬ 
ated with the creep of the piezoelectric scanner. By con¬ 
trast, u y vg (x) exhibits a sharp kink at the TB [Fig. 6(e)]. 
These features are consistent with the model shown in 
Fig. 6(a). It should be noted that du y vg (x)/dx shown in 
Fig. 6(f) is almost completely constant in both domains, 
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FIG. 6 : (a) A schematic top view of the atomic arrangement near the TB of FeSe (not in scale). Note that an atomic 

periodicity along the x direction is hardly affected by the TB. (b) A constant-current STM image taken over a field of view of 
180 nmx90 nm on a grid of 4096x2048 pixels. The set-up conditions for imaging were V s = +95 mV and It = 100 pA. Inset: 
Fourier-transformed STM image taken in the left domain of the main panel. A peak at q x is sharp and well isolated from other 
features, guaranteeing that q x ■ u(r ) can be treated as a spatially varying phase of the q^-modulations. (The same is true for 
q y .) (c) The x component of u(r ) averaged over the y axis. Thick solid line and thin dashed line denote the data taken by the 
forward (left to right) and backward (right to left) scans, respectively. The symmetric hysteretic behavior between the forward 
and backward scans indicates that u™ s (x) is governed by the creep of the scanner. No anomaly is observed at the TB. (d) 
The x derivative of u“ vg ( x). Spike-like features are associated with the point defects in the image, (e) The y component of 
u(r ) averaged over the y axis. Thick solid line and thin dashed line denote the data taken by the forward (left to right) and 
backward (right to left) scans, respectively. Since the y direction is the slow-scan direction, the effect of the creep is small. A 
clear kink is observed at the TB. (f) The x derivative of Uy VS (x). 


indicating that the TB-induced strain to the lattice is 
negligibly small. 

The observed value of dUy VS (x)/dx ~ —1.1 x 10 “ 2 
in the right domain means that the angle /? defined in 
Fig. 6 (a) is +0.63°. This means that orthorhombic dis¬ 
tortion (b — a)/(b+ a) « 2.8 x 10 ~ 3 , being consistent 
with the X-ray result (33j. Even if there were an addi¬ 
tional lattice distortion associated with the TB, it should 
be much smaller than this tiny orthorhombic distortion 
which we have clearly detected. 


Appendix B: Asymptotic forms of the order 
parameter derived by the Ginzburg-Landau theory 


We derive asymptotic forms of the order parameter 
far from TBs using the Ginzburg-Landau (GL) theory. 
We consider the GL free-energy functional for tetragonal 
symmetric systems [30] as an expansion in the isotropic 
s-wave component Aj so and the four-fold d- wave compo¬ 
nent A 40 of the order parameter: 


F GL [A iso ,A 40 ] = Jdvl Y, [ay(T)\A^\ 2 + b^\AX + K^\VA^\ 2 ]+^\A iso \ 2 \A^ 


/i,=iso,40 


+? (A* 2 A 2 0 + A 2 so A H) + ^ p a A iso )*(d a A 40 ) - ( d b A iso )*(d b A H ) + c.c.] \ , (Bl) 


K 






























10 


where we have neglected the vector potential as it does 
not play an important role in our discussion. The co¬ 
efficients b/j,, K^, and K are positive and a M (T) = 
a M (T/T CjU — 1) with positive a The differential oper¬ 
ator V = (d a ,db) is defined according to the crystal axes 
a and b. As in Ref. [30|, we assume 72 > 0, so that the 
free energy is minimized at ip = ±7r/2, and the time- 
reversal-symmetry-broken s ± id state is stabilized when 
both Aj so and A 4? j are finite. 

The effect of orthorhombic distortion is taken into ac¬ 
count by adding the following term to the free-energy 


functional [30| : 

F e = ce j dV(A* so A^ + A iso A: 0 ), (B2) 

where c is a positive parameter and e = e aa — £bb is the 
parameter of the orthorhombic lattice distortion. The 
total free energy for a uniform state in the bulk is then 
given by 


Fch + Fe = J2 (a,\A,\ 2 +b,\AX) 

fi—iso,4(f) 

+ 7i|Ai so | 2 |A 40 | 2 + 7 2 |A iso | 2 |A 40 | 2 cos( 2 ^) + 2 ce|A iso ||A 40 | cos^, (B 3 ) 

where A M = \A li \e lv i i and the relative phase ip = <p 4 ^ — ipi so - If c|e| > 272 1Aj so 11A 4 ^|, then the free energy is minimized 
at ip = 0 for e < 0 and at ip = tt for e > 0. In the following discussion we assume that this inequality is satisfied and 
the time reversal symmetric s ± d state is realized in the bulk. 

Next, we consider a TB located at x = Xq along the y axis, where the x and y axes are rotated by 45° from the 
crystalline axes, x = (a — b)/V2 and y = (a + 6 )/v2. The orthorhombic lattice distortion parameter e changes its 
sign across the TB. We assume e(x) =F |e| for x —> ±00, so that the s ± d state is realized in x — ¥ ±00. Near the TB 
where e(x) is small, the s ± id state is favored. Then, the area density of the total free energy is given by 


/gl + ft = / dx 


(S/i | A m 1 2 + A^| 4 + K^ldxA^ 2 ) + 71I Ai so | 2 | A 4? i|^ 

. fi= iso,40 

+ ■y (^i*so^40 + ^?so^4 2 ) + Ce ( a; )(^i*so^40 + AisoAl^) 


(B4) 


Let us first assume that Ai SO is a real and uniform order parameter while A 4 ^ changes its sign across the TB, as 
shown in Fig. 4(a). If we restrict A 4 ^ to be real, then A 4 ^ varies over the coherence length [3d] 


£ 


K^cj, 


■ 66 40 |A ^|2 + (7i+72) | A . 


G&4</> 


(B5) 


where | A^ lk | is the amplitude of A 4 ^ in the bulk. However, we expect that time-reversal symmetry should be locally 
broken at the TB. Thus, we allow A 4 ^ to be complex, A^x) = \ A^{x)\e ltf, ^ x \ With this order parameter, the total 
free energy is given by 


/gl + ft = dx 


®iso|Aiso| 2 + I>iso|Aj SO | 4 + a 4 0 1A 4 0 1 2 + 1> 4 01A 4 01 4 + 7:L|A iso | 2 |A 4 4,| 2 

+ 7 2 |A iso | 2 |A 4 0| 2 cos(2 ip) + 2ce(a; - xo)|A isc ,||A 4 ^| cos 73 
+ K 44> [(5 x |A 4 ^|) 2 + |A 4 0| 2 {d x (p) 2 


I- 

/gl + ft is 
K, 


(B6) 


= - | A buik| [ 72 |A iso ||A^ lk | sin(2 p) - c|e| sin<^] . 
1^40 I 

(B7) 


In the bulk region (x — Xo £) where | A 4 ^| = 0 and 
e(a;) = — |e|, the GL differential equation to minimize 
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Since p <C 1 far away from the TB, we can linearize the 
differential equation and find the relative phase to decay 
as ip oc exp (— x/£) with the characteristic length 


£ = 




\ |A iso |(c| e |- 2 72 |A iso ||A^k|)- 


(B8) 


The characteristic length diverges when approaching the 
phase boundary, where c|e| = 2 ^ 2 1 Ai SO 11 A 4 ^ |, between 
the time reversal symmetric sic? state and the time- 
reversal-symmetry-broken s ± id state. 


Finally, we consider double TBs at x = i|xo|, where 
£ -C |xo| < £. We assume e > 0 between the double TBs 
and e < 0 otherwise. At the center x = 0 between the 
double TBs, we can set | A 40 I = |A^ lk | because |xo| £. 
With this approximation, the GL differential equation to 
minimize /ql i ft for |x| -C |xo| is 

K^dlp = [ 7 2 |A iso ||A^ lk | sin(2^) + c|e| sin^] . 

(B9) 


Integration of the differential equation yields 


K^{d x p) 2 


l ^ui| I [72 1 Also 11 A 4 ^ lk | cos(2 ip) i 2c|e| cos ip- 7 2 |A iso ||A^ lk | cos(2 pff) - 2c|e| cos<p 0 ], 

l A 4 0 I 


(BIO) 


where the integration constant is determined from the 
conditions d x p{x = 0) = 0 and p[x = 0) = ipg. Since 
we assume the distance between the TBs is in the range 
|xo| < £, the relative phase does not reach n at x = 0 , 
i.e., ipo < n. For po — p -C 1 near x = 0, the differential 
equation (IB 101) has the solution 


<p{x) =tp 0 - 



(Bll) 


For the model order parameter shown in Fig. 5(d), we 
have determined po and £0 by the continuity condi¬ 
tion at x = ±|xo|/2, that is, by imposing that p = 
7 r — (7r/6)sech(x/£) and p(x) in Eq. (IB11I) are smoothly 
connected at x = ±|x 0 |/2. We note that different choices 
of the connecting position yield little change in the value 
of p. 
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